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Abstract. The paper introduces a general framework for derivation of continuum equations governing meso-scale dynamics 
of large particle systems. The balance equations for spatial averages such as density, linear momentum, and energy were 
previously derived by a number of authors. These equations are not in closed form because the stress and the heat flux cannot 
be evaluated without the knowledge of particle positions and velocities. We propose a closure method for approximating fluxes 
in terms of other meso-scale averages. The main idea is to rewrite the non-linear averages as linear convolutions that relate 
micro- and meso-scale dynamical functions. The convolutions can be approximately inverted using regularization methods 
developed for solving ill-posed problems. This yields closed form constitutive equations that can be evaluated without solving 
the underlying ODEs. We test the method numerically on Fermi-Pasta-Ulam chains with two different potentials: the classical 
Lennard- Jones, and the purely repulsive potential used in granular materials modeling. The initial conditions incorporate 
velocity fluctuations on scales that are smaller than the size of the averaging window. The results show very good agreement 
between the exact stress and its closed form approximation. 
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1. Introduction. Particle systems governed by Newton's ordinary differential equations (ODEs) are 
common in physics, engineering and computational biology. Such systems could represent either classical 
molecular models, or discretizations of continuum mechanical partial differential equations (PDEs). These 
ODE systems are difficult to simulate directly because of their large size and stiffness. Popular explicit 
solvers such as Verlet method require small time steps, which places severe restrictions on the length of 
the simulated time interval. This necessitates the development of new methods for reducing computational 
complexity. 

Often, the ODE solutions are of secondary interest compared with various space-time averages. Examples 
of the latter are density, velocity, stress, deformation gradient, and energy. The averages can be always 
simulated directly, but a better option is to formulate an approximate, continuum mechanical type model 
that describes the evolution of the averages. Such a model can be simulated at a lower cost than the 
underlying ODE system. 

The theory of space-time averaging for particle systems was developed by several authors, starting with 
Irving and Kirkwood '8] and Noll [19]. Several decades later. Hardy [7] and Murdoch and Bedeaux [16], 
[17j . jl8j . [TS] developed the theory further and derived the governing balance equations. The averaging in 
these works is done as follows. First, one selects the primary averages that would describe the meso-scopic 
state of the ODE system. These may be, for example, average density, velocity, deformation map, kinetic 
energy etc. Next, differentiating in time and using the ODEs, one obtains the governing balance equations 
that model meso-scale behavior of the particle system. The fluxes (or secondary averages) in these equations 
are given by explicit functionals of the ODE solutions. This establishes a connection between the fine scale 
ODE model, and the meso-scale PDE model. 

While important for clarifying the relationship between micro- and meso-scale phenomena, results of this 
type do not provide a continuum model in the true sense of the word, because one still has to solve the ODEs 
to evaluate the fluxes. Therefore, the balance equations in |TB] are not in closed form. In classical continuum 
mechanics, the constitutive equations express stress and heat flux in terms of velocity gradient, deformation 
gradient, and temperature. In the above referenced theories, the average deformation and temperature are 
not sufficient for evaluating fluxes, because one still has to recover the positions and velocities of all particles. 
As a result, the complexity of the meso-scale models in [T^ and [7] is about the same as the complexity of 
the ODE model. 
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In this paper, we propose a method for approximating fluxes in terms of the primary averages. The 
approximations have the same low complexity as the classical constitutive equations, but unlike these equa- 
tions, our approximations combine explicit formulas and numerical algorithmic prescriptions. From our point 
of view, any constitutive equation will be ultimately realized as a computational method. The accuracy and 
efhciency of this method are the main factors that determine the quality of the constitutive approximation. 
Therefore, instead of trying to construct short explicit equations, we sought a computational scheme that 
(i) does not require solving the ODE system; (ii) clearly and consistently incorporates the micro-scale force 
equations; and (iii) reproduces the exact fluxes with reasonable accuracy. 

The stress in our constitutive equations depends on density and momentum in a non-local and non- 
linear manner. The non-locality makes our theory similar to the peri-dynamical formulation of continuum 
mechanics |22j . The difference between our approach and phenomenological peri-dynamics is that our method 
clearly links the micro-scale features of the particle dynamics and meso-scale constitutive equations. On the 
other hand, our work differs from the recent paper |13j where the peri-dynamic stress and energy flux are 
given as exact functions of particle positions and velocities. In our theory, the stress can be evaluated without 
solving the ODE system. 

The main idea behind our approach is as follows. The primary averages are essentially non-linear integral 
operators acting on particle positions and velocities. These operators can be written as (linear) convolutions 
of the "window function" and certain functions of velocities and positions. In this way, each primary average 
is related to a function of the micro-scale dynamics. For example, density corresponds to the Jacobian of 
the inverse Lagrangian deformation map, and linear momentum corresponds to the product of this Jacobian 
and micro-scale velocity. The convolution operator is typically invertible, and thus micro-scale quantities 
can be, in principle, recovered from the averages. However, the deconvolution problem is unstable (ill-posed) 
so that small perturbations of the averages can produce large perturbations in the recovered functions. 111- 
possed problems are well studied in the literature (see e.g. [3 [9l O [TH [21] ) in both continuous and discrete 
settings. A discrete version of the convolution integral equation is an ill-conditioned linear system. Such 
systems and related rank-deficient systems are treated in detail in the book [B]. The strategy for solving 
ill-posed problems is to approximate the exact problem by a well-posed regularized problem depending 
upon a parameter. By varying this regularization parameter one obtains more accurate but less stable 
approximations. Loosely speaking, the effect of regularization is to smooth the exact solution and filter out 
higher frequency oscillations. The amount of smoothing and filtering depends of the choice of the method 
and the value of the regularization parameter, but some details of the exact solution are always lost. Despite 
this, it is always possible to produce a regularized deconvolution that is much closer to the exact micro-scale 
function than the corresponding average. The difference between the average and the regularized solution 
can be quite dramatic. The convolution, whose kernel is associated with the mesoscopic length scale, smears 
smaller details to such an extent that it is often impossible to recognize them by inspecting the graph of the 
average. A well chosen regularization performs a triage of length scales below the mesoscopic length: the 
smallest are filtered out, and the rest are recovered. 

Let us now describe the main steps of the method. 

(i) The starting point is a fine scale model: an ODE system of Newton's equations. We limit ourselves 
to the case of pairwise, short-range interaction forces that may depend on the relative positions and 
velocities. The system size N is very large. A typical interparticle distance is characterized by a 
small parameter 

e = 7V-i/^ (1.1) 

where d is the physical space dimension (usually 1, 2, or 3). The particle masses and forces are 
scaled by e. The purpose of the scaling is to satisfy the following natural requirements. As iV oo, 
the total mass of the system should remain fixed, and the total particle energy should be either 
fixed, or at least bounded independent of N. 

The scaled ODE systems with increasing N form a family of discrete models representing a single 
continuum system at different levels of micro-scale resolution. To prescribe initial conditions, we 
first fix the initial velocity interpolant. Then, given N, the initial particle velocities are generated 
by discretizing the interpolant on the uniform mesh of size sL. The initial positions are the nodes 
of the same mesh. 
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(ii) Choose spatial mesoscale resolution parameter rj. Fix a function ip with integral equal to one and 
compact support (non-compactly supported functions such as Gaussian are also possible). Then 
scale this function by a factor of rj and define the window function 



The window function is used to set up meso-scale averaging, 
(iii) Choose the primary mesoscopic averages (e.g. density, velocity, internal energy density, temper- 
ature). The collection of the primary averages represents the state of the particle system at the 



(iv) Write down the exact balance equations for the primary variables. Determine which fluxes in these 
equations require closure, and list all the microsopic quantities (e.g. interpolants of particle positions 
and velocities, Jacobians) that need to be reconstructed to achieve closure. 

(v) Verify that the every required micro-scale quantity can be (approximately) reconstructed from the 
chosen primary variables. In each flux, replace the exact microscopic quantities with their regularized 
deconvolution approximations, thereby obtaining approximations of the microscale quantities in 
terms of the primary averages. 

We tested the method numerically on one-dimensional Hamiltonian systems with short-range pair po- 
tentials. In the physics literature such systems are known as Fermi-Pasta-Ulam (FPU) chains. We consider 
two potentials: the classical Lennard- Jones, and another potential similar to the Hertz potential of granular 
acoustics. Assuming that the meso-scale state of a system can be described by the density and linear mo- 
mentum, we provide the corresponding balance equations and derive constitutive equations for the stress. 
The exact stress and exact primary averages are produced by the direct simulations with 10,000 particles. 
The approximate stresses are obtained by using deconvolution and then substituting into the formulas for 
the exact stress. Then we compare the exact and approximate stresses rendered on the meso-scale mesh 
with 500 nodes. The results show that the approximation agrees very well with the exact stress. 

The present article extends and improves the method introduced in |20j . In |23| , some of the tools from 
PO] are applied to the discrete models of fluids. In both papers, the suggested deconvolution algorithm 
was the classical Landweber iteration [1], [lOj . Increasing the number of iterations n increases accuracy 
but generally decreases stability. The zero-order approximation (n = 0) consists of replacing micro-scale 
quantities with their averages. This zero-order closure was studied in detail in |20j . In |23j we also used the 
first- and second-order approximations. Numerical experiments show that low-order closures work well when 
the initial velocity has small fluctuations, and the dynamics is nearly isothermal, meaning that the energy 
of velocity fluctuations is much smaller than the potential energy. 

The Landweber iteration is simple and useful for modeling, but has a slow convergence rate (see [S]). For 
initial conditions with high fluctuations, a large number of iterations may be needed to achieve a reasonable 
accuracy. In this work we use different techniques: regularization by discretization [5] and truncated singular 
value decomposition (SVD) (see, e.g. [B]). Both methods are non-iterative. Regularization by discretization 
is straightforward: the integral is approximated by a numerical quadrature, and this eliminates accumulation 
of the spectrum to zero. In the truncated SVD method, the exact solution is represented in the basis of 
singular vectors. The regularized approximation is generated by discarding the components corresponding 
to the smallest singular values. Using SVD yields additional computational savings, since the convolution 
kernel is dynamics- independent. The SVD of the kernel can be pre-computed and used repeatedly with 
different ODE systems. 

Recently, a deconvolution approach was used in large eddy simulation (LES) of turbulence [I], [1], [H], 
[TT] . In these works, deconvolution was used to approximate quadratic functions of velocity fluctuations by 
an operator acting on the average velocity. The present work differs from LES in several respects. The flrst 
difference is in the structure of the averaging operators. In LES, the average velocity depends linearly on the 
micro-scale velocity, while in our work this dependence is non-linear. This non-linearity makes it possible to 
handle the general ODE flows with non-constant Jacobians. The second difference is in the modeling. We 
provide the connection between Newtonian particle mechanics on the one hand, and continuum theories on 
the other hand. In LES, the objective is to simulate large scale features of flows governed by Navier-Stokes 
equations. The third difference is that the papers on LES do not systematically address ill-posedness, and 




(1.2) 



mesoscale. 
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do not make much use of the available results on ill-posed and inverse problems. The Gaussian kernel, one 
of the most popular in LES, is not the best in terms of stability of reconstruction, because the degree of 
ill-posedness depends on the smoothness of the kernel in the Sobolev scale: the smoother the kernel, the 
more unstable the reconstruction problem. For this reason, we use piecewise polynomial continuous kernels, 
which leads to a mildly ill-posed problem. 

The paper is organized as follows. In Section 2 we describe a general multi-dimensional microscopic model 
and provide the exact balance equations for the averages, following [IS], [IS]. In Section 3 we develop general 
multi-dimensional integral approximations of averages, and describe the use of regularization for approximate 
deconvolution. This section is the central section of the paper. Section 4 contains the formulation of the scaled 
ODE equations of FPU chains. In Section 5 we derive closed form balance equations of mass and momentum 
for such chains and provide the constitutive equations. Section 6 contains the results of computational tests. 
Finally, conclusions are given in Section 7. 

2. Microscale equations and mesoscale spatial averages. 

2.1. Scaled ODE problems. We work with classical Newton equations of point particle dynamics. 
The same equations may arise as discretization of the momentum balance equation for continuum systems. 
Consider a system containing ^ 1 identical particles, denoted by Pi. The mass of each particle is 
where M is the total mass of the system. Suppose that during the observation time T, Pi remain inside a 
bounded domain in W^^ where d is the physical space dimension, usually 1, 2 or 3. The positions <7j(i) and 
velocities Vi{t) of particles satisfy a system of ODEs 

= (2.1) 

^v, = f. + ft'\ (2.2) 



subject to the initial conditions 



q,(0) = a;„ v,{0) = v1. (2.3) 



Here f'^^^*^ denotes external forces, such as gravity and confining forces. The interparticle forces = /j^ , 
where fi^ are pair interaction forces which depend on the relative positions and velocities of the respective 
particles. 

We are interested in investigating asymptotic behavior of the system as TV — )■ oo. Thus it is convenient 



to introduce a small parameter e by (1.1) characterizing a typical distance between neighboring particles. 



As e approaches zero, the number of particles goes to infinity, and the distances between neighbors shrink. 



Consequently, the forces in (2.2) should be properly scaled. The guiding principle for scaling is to make the 
energy of the system bounded independent of A'^, as iV — cxd. In addition, the energy of the initial conditions 
should be bounded uniformly in N . 

As an example of scaling, consider forces generated by a finite range pair potential U{S^)^ where ^ is the 
distance between the particles. Suppose that each particle interacts with no more than a fixed number of 
neighbors. This implies that there are about N interacting pairs. If the system is sufficiently dense, and 
variations of particle concentrations are not large, then a typical distance between interacting particles is on 
the order N'^^'^L = eL. The resulting scaling 



makes the potential energy of an isolated system bounded independent of N. Kinetic energy will be under 
control provided the total energy of the initial conditions is bounded independent of N. If exterior forces 
are present, they should be scaled as well. 



Remark. Superficially, the system ( 2.1 ), ( 2.2 ) looks similar to the parameter-dependent ODE systems studied 
in numerous works on ODE time homogenization (see e.g. [21j and references therein). In the problem 
under study, e depends on the system dimension N, while in the works on time-homogenization and ODE 
perturbation theory, the system size is usually fixed as e — > 0. 



2.2. Length scales. We introduce the following length scales: 

- macroscopic length scale L = (iiam(f2); 

- microscopic length scale eL\ 

- mesoscopic length scale 77L, 

where < ry < 1 is a parameter that characterizes spatial mesoscale resolution. This parameter is chosen 
based on the desired accuracy, the computational cost requirements, and prior information about initial 
conditions and ODE trajectories. 

The computational domain f2 is subdivided into mesoscopic cubic cells C^, 13 = 1,2, ... ,i? with the 
side length on the order of rjL. The centers xp oi Cp are the nodes of the meso-mesh. The number of 
unknowns in the mesoscopic system will be on the order of B. For computational efficiency, one should have 
B <ti N. This does not mean that 77 is close to one. In fact, it makes sense to keep rj as small as possible in 
order to have an additional asymptotic control over the system behavior. Decreasing rj will in general make 
computations more expensive. 

2.3. Averages and their evolution. To define averages we first select a fast decreasing window 
function satisfying J ip{x)dx = 1. There are many possible choices of the window function. We assume 
that -0 is a compactly supported, continuous, differentiable almost everywhere on the interior of its support, 



and non-negative. Next, define the window funciton ijjjj by (1.2 1. 

Once the window function is chosen, one can generate averages of micro-scale dynamical functions (|16j. 
|15j). The mesoscopic average density and momentum are given, respectively, by 

M ^ 

r{t,^) = J;^T.M^-1^(t)), (2.5) 

i=l 

p^v^t, ^) - f E I'M- (2-6) 

The meaning of the above definitions becomes clear if one considers ip = {cd)^^x{^)j where x is a charac- 
teristic function of the unit ball in M'^, and Cd is the volume of the unit ball. Then 

CdTj"- JM ^-^ \ rj 

The sum in the right hand side gives the number of particles located within distance 77 of a; at time t. 
Multiplying by M/N we get the total mass of these particles, and dividing by Cdi]''' (the volume of ry-ball) 
gives the usual particle density. 



Differentiating (2.5), (2.6| in t, and using the ODEs (2.1|, (2.2) one can obtain [16] exact mesoscopic 
balance equations for the primary variables. For example, for an isolated system with {f^^^^^ — 0), mass 
conservation and momentum balance equations take the form: 

^tp" -I- div(p''tJ'') = 0, (2.7) 

dt{]fv'') + div (p^U" ® U") - div T" = 0. (2.8) 
The stress T"* = Tj'^^ + T|'^„^j [H], where 

T^l^^it, x) = -Y, - v'^it, X, )) {v, - v'\x, t))ij,^{x - q,) (2.9) 

is the convective stress, and 

r''(t, a;)(„t) = V /,y (g) {Qj - Qi) [ i'r, {s{x - Qj) + (1 - s){x - q^)) ds 



(2.10) 



is the interaction stress. The summation in (2.101 is over all pairs of particles («, j) that interact with each 
other. 

Discretizing balance equations on the mesoscopic mesh yields a discrete system of equations, called 
the meso-system, written for mesh values of p^^, ip^v'^)^ and Tj^. The dimension of the meso-system is 
much smaller than the dimension of the original ODE problem. However, at this stage we still have no 



computational savings, since the meso-system is not closed. This means that mesoscopic fluxes such as (2.91, 



(2.10) are expressed as functions of the microscopic positions and velocities. To find these positions and 



velocities, one has to solve the original microscale system (2.11, (2.2|. To achieve computational savings we 
need to replace exact fiuxes with approximations that involve only mesoscale quantities. We refer to the 
procedure of generating such approximations as a closure method. This closure-based approach has much 
in common with continuum mechanics. The important difference is that the focus is on computing, rather 
than continuum mechanical style modeling of constitutive equations. 

3. Closure via regularized deconvolutions. 

3.1. Outline. Our approach is based on a simple idea: the integral approximations of primary averages 
(such as density and velocity) are related to the corresponding microscopic quantities via convolution with 
the kernel ipri. Therefore, given primary variables we can (approximately) recover the microscopic positions 
and velocities by numerically inverting convolution operators. The results are inserted into equations for 
secondary averages (or fluxes), such as stress in the momentum balance. This yields closed form balance 
equations that can be simulated efficiently on the mesoscopic mesh. 

3.2. Integral approximation of discrete averages. To exploit the special structure of primary 
averages, it is convenient to approximate sums such as 

_ 1 ^ 1 \n\ ^ 



by integrals. The sum in (3.1 1 resembles a Riemann sum for \ri\ g'ip,^{x — •), where VL is partitioned into N 



cells of volume with one particle located inside of each cell. However, because of the motion of particles. 



(3.1) is not in general a Riemann sum. Indeed, to interpret this sum correctly, one must exhibit a partition 
of into cells of equal volume where each cell contains exactly one particle. Such a partition may not exist. 
Indeed, in one dimension, the domain is an interval, say (0, L) and the cells are intervals of length L/N. 
Thus, if the closest neighbors of a given particles are less than L/N apart, than the desired partition does 
not exist. For two- and three-dimensional domains, it may be possible to use more general partitions, but 
for particles that are spaced non-uniformly, the shapes of these cells may be quite far from slightly deformed 
rectangles, which would make it difficult to estimate the accuracy of the resulting integral approximation. 

A more systematic way to generate integral approximations is to use a microscopic flow map interpolant 
and the associated Jacobian describing local volume changes. Let q{t, X),v{t,q) be suitable position and 



velocity interpolants, associated with the system (2.1), (2.2). At < = these interpolants satisfy 

q{0,X,) = q°, i3(0,9(0,X,)) = t;^", 

where Xj, j = 1,2, . . . , N are points of e-periodic rectangular lattice in Q. At other times, 

q{t,X^) = qJt), v{t,q{t,X,))=v,{t). 



Then we can rewrite (3.1 ) as 

N 



1 101 



I I j=i 



where \Q\ denotes the volume (Lebesgue measure) of fl. Eq. (3.2 ) is a Riemann sum generated by partitioning 
O into N cells of volume |0|/iV centered at Xj. This yields 



l^Jj {v{t, q{t, X)),qit, X)) ij^ix - qit, X))dX, (3.3) 



up to discretization error. Now suppose that the map q{-, X) is invertible for each i, that is X = q ^{t,q). 
Changing the variables in the integral y — q{t, X) we obtain a generic integral approximation 



5'' = j^g {v{t,y).y)'<IJr,{x - y)J{t,y) dy, 



(3.4) 



where 

J=|detVq"^|, (3.5) 

up to discretization error. For reader's convenience, we list the integral approximations of the average density 
and momentum. 

M ^ 

p^{t,x)^-Y.^,{x-qM (3.6) 



N 



M 

Pi 

M 

M 



/ ^^{x-~q{t,X))dX 

/ tl}r,{x~y)J{t,y)dy. 
Jn 



M ^ 

P^v^it, x)^-J2 v^(t)M^ - qm (3.7) 



N 
M 

¥\ 

M 

W\ 



I v{trq{t,X))i>,,{x~~q{t,X))dX 
ipri{x - y)v{t,y)J{t,y)dy. 



Note the linear convolution structure of the y-integrals in (3.6 1 and (3.7 1. It is also worth noting that these 
equalities are exact if the interpolants are piecewise linear. In that case, the discrete sums are exact integral 
quadratures. 

3.3. Regularized deconvolutions. 

3.3.1. General considerations. Define an operator Rjj by 



RrjifK^) = y ^v{x~ y)f{y)dy. 



To simplify exposition, suppose that i?,, is injective. In that case, there exists the single-valued inverse 
operator that we call the deconvolution operator. Since Rr^ is compact in L^(n), the inverse operator 

is unbounded. Therefore, small perturbations of the right hand side can lead to large perturbations in the 
computed solution. Reconstructing / from the knowledge of -Rr;[/]) is a classical example of an unstable 
ill-posed problem. Such problems are well investigated both analytically and numerically (see, e. g. O 
[9l [HI 1211 13] ) . Many solution techniques are currently available: Tikhonov regularization, iterative methods, 
reproducing kernel methods, the maximum entropy method, the dynamical system approach and others. In 
the sequel we use the notation for a regularized approximation of the exact inverse operator. 

Recently, the classical Landweber iterative deconvolution method |3], [TO] has attracted attention as a 
means to achieve sub-filter scale resolution in large eddy simulation of turbulence [I], [1]. In the simplest 
version of the method, approximations g„ to the solution of the operator equation 



are generated by the formula 



i?„[.g]=r' (3.8) 



gn = Y.{I-R^r9\ 50=5". (3.9) 

fc=0 
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The number n of iterations plays the role of regularization parameter. In (3.91, / denotes the identity 
operator. 



Another classical method is Tikhonov regularization |24j . where the solution of (3.8 1 is approximated by 
ga that solves 

RM+aC[g^]^g\ (3.10) 

Here a > is a regularization parameter, and C is a stabilizing operator. In practice, C can be an identity, 
or a a suitable differential operator such as Laplacian. 

3.3.2. Computational implementation of deconvolution. The discrete version of the integral 



equation (3.81 is a linear system 

Ag = g. (3.11) 

The solution g is a discrete micro-scale quantity to be reconstructed, and g is the discretization of the 
corresponding average. To achieve computational efficiency, it is natural to resolve the average (a meso-scale 
quantity) on a coarse mesh with the size tied to the meso-scale. The solution (a microscopic quantity) could 
be rendered on the fine mesh with size eL. This choice of meshes seems to be the most natural for balancing 
cost and accuracy. 

Other mesh combinations can be chosen as well. The least expensive option is to coarsen the discretiza- 



tion of g and solve (3.11 1 entirely on a coarse mesh. In that case, the matrix A is square, and its dimension is 
determined by the chosen number of the coarse mesh points. The operation count of deconvolution becomes 
independent of TV, but this efficiency comes at the price of introducing too much error. Numerical experi- 
ments produced significant artifacts, so that the meso-scale details of g could not be reliably reconstructed. 

Another possibility is to use a fine mesh for both g and g. In that case the average would have to be 
interpolated from the coarse to the fine mesh. The computational cost of deconvolution scales as 0{N'^), 
but this dos not necessarily improve the reconstruction quality, since using a finer discretization increases 
the condition number of A. 

Ultimately, we chose the more natural two-mesh discretization, whereby g is an iV-vector, g is a D- 
vector, with B < D N . Recall that B is the number of the mesh nodes associated with the mesh size rjL. 
Choosing a finer coarse mesh with D nodes enabled us to better resolve the details of size smaller than rjL. 
Such details may be completely smeared by averaging. In that sense, the length scale associated with D can 
be called a sub-filter scale. 

Working with two different meshes, as opposed to using the same mesh, introduces several difficulties. 



The matrix A in that case is rectangular, and the system (3.11) is under-determined. Therefore, (3.111 has 
many solutions, even when the original integral equation (3.8) is uniquely solvable. The general solution is 
a sum of the particular solution g^ orthogonal to the null space of A and an arbitrary vector from the null 
space. In the absence of a priori information on the structure of the null space, it is natural to use g^. Thus 
we set 

g+^A^z, 

where 2; is a Z?-vector to be determined. Assuming that A has full rank and AA^ is invertible (invertibility 
depends only on the choice of the window function tp and can be verified prior to running any dynamic 
simulations), we can set z — {AA'^)~^g, and 



g+ = A' (AA' )-'g. (3.12) 



It is easy to check that g+ is orthogonal to the null space of A. In (3.12), {AA^)~^ denotes either the 
exact inverse, or a suitable regularized approximation. Typically, singular vectors associated with smaller 
singular values of A oscillate with higher frequency than the vectors associated with larger singular values. If 
this holds, then the solution component along the null-space of A is highly oscillatory, while the component 
orthogonal to the null space is relatively smooth. Therefore, by using g^ and not some other solution, we 
incorporate additional filtering. This can be useful for taming noise. 



An example of deconvolution is shown in Figs. 3.1 and 3.2 In the first of these figures we show the exact 
solution g. It is constructed by first choosing a profile to be reconstructed (left panel), and then adding 
noise (right panel). The graph in the left panel contains a meso-scale feature (trapezoidal impulse in the 
center), and a sub-filter scale feature (smaller triangular impulse on top of the trapezoid). The noise contains 
multiple features on smaller length scales. 

The right hand side vector g, computed by applying ^ to g in the right panel of Fig. |3.1[ is shown 



in Fig. 3.2 



All features except the largest appear to be smeared. The right panel in Fig. |3.2| shows the 

The noise is largely filtered out but the sub-filter feature is clearly 



reconstruction computed using (3.12) 
visible. 





Fig. 3.1. Left panel: meso-scale and sub-filter scale features; right panel: exact solution with a uniformly distributed noise 
added 





Fig. 3.2. Left panel: the average (right hand side of the integral equation); right panel: reconstructed approximate solution 



Many regularization methods (Tikhonov, Landweber, truncated SVD) can be conveniently written in 
terms of SVD and spectral filter functions (see e.g. [9J). This is useful for both discrete and continuous 
ill-posed problems fB^. Since both i?,, and A are independent of the microscopic dynamics, the SVD can be 
pre-computcd and used with different ODE systems. 

Suppose that rank(y4) = D. Let Cj, j — 1, . . . , 13 denote the non-zero singular values of A, and G R^, 

e be the corresponding singular vectors. By the standard properties of SVD, 



(3.13) 



Because the continuum problem (3.8) is ill-posed, the singular values are spaced without gaps, and the 



condition number of A is large. Recall that condition number can be expressed as the ratio of the largest 
and smallest singular values. In our case, the largest singular value of A is close to one. Therefore, the 
condition number is approximately equal to the reciprocal of the smallest singular value. The choice of the 
solution method depends on the condition number and the relative level of noise in g and A. The guiding 



principle it to produce an approximation that would be as close as possible to (3.121, without incurring 
instability. Discretization itself is a mild form of regularization of the original integral equation. Indeed, 
in the discretized problem, the smallest non-zero singular value is a finite distance away from zero, while 
in the continuum case zero is an accumulation point of the spectrum. Consequently, if the error in g and 

9 



the condition number of A are small enough, one can use (3.12) with no additional regularization. In our 



case, the integral approximations of the averages are exact when the interpolants are suitably chosen. The 
only numerical error in the right hand side is the round-off error. Therefore, exact inversion will work when 
the condition number of A is much smaller than the reciprocal of the machine precision. In practice this 
means that condition numbers smaller than about 10^ can be safely handled in this way. For larger condition 



numbers, exact inversion in (3.121 may have to be replaced by a suitable regularized approximation. 



First we describe the SVD-based implementation of the exact solution formula (3.121. Write 



D 



D 



9 = II-9j^j' 



(3.14) 



where the coefficients Zj have to be determined. To obtain the last equality, (3.13) was used. Substituting 
(3.14) into (3.11) and using orthogonality of we deduce 
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(7i 



(3.15) 



which represents (3.12) in the basis consisting of singular vectors. 



As was explained earlier, formula (3.15) may be used when condition number of A is much smaller than 



the reciprocal of the noise level in the right hand side. Therefore, (3.15) becomes unstable when the SVD 



contains singular values comparable to the machine precision. In that case, as an additional regularization 
we use the truncated SVD 6 . In this method, the components corresponding to the smallest singular values 
are discarded. The regularized solution is computed by the formula 



where the filter function (j) is defined as follows. 

(/-(cTj) = 



D 

E 



1 if a J > a* 
if (T,- < a* 



(3.16) 



(3.17) 



In the above equation a* is a cut-off value (equal to the machine precision in the present case). 

4. FPU chain equations. In this section, the general method outlined above is detailed in the case 
of one-dimensional Hamiltonian chain of oscillators that consists of N identical particles. The domain 17 is 
an interval (0, L). Particle positions, denoted by qj = qj{t), j — 1, . . . , N , satisfy 

< qi < q2 < . ■ ■ < qN < L 

at all times, i.e. the particles cannot occupy the same position or jump over each other. Next, define a small 
parameter 

1 



and microscale step size 



The interparticle forces 



fjk 



1] - Ik jj, fkj -<ik\ 



(4.1) 



(4.2) 



are defined by a finite range potential U. 



10 



Each particle has mass m — M/N ~ Me, where M is the total mass of the system. Particles have 
velocities denoted hy Vj, j = 1, . . . , iV. Writing the second Newton's law as a system of first order equations 
yields the scaled microscale ODE system 



eMv,=fj, j = l,...,7V 



(4.3) 



subject to the initial conditions 



(4.4) 



5. Integral approximation of stresses for particle chains. Mesoscopic continuum equations. 



In the one-dimensional case stress is a scalar quantity, and (2.9), (2.101 reduce to, respectively, 

^ M 



(5.1) 



and 



.7 = 1 ^0 



(1 — s)qj)ds 



(5.2) 



The sum in (5.2| is simplified compared to the general expression (2.10), since we have exactly iV — 1 
interacting pairs of particles. 



To obtain integral approximations of stresses, we define interpolants q,v, as in Sect. 3.2 Repeating the 
calculations we get 



M 



{v{t,y) - v''{t,x)) ipr,ix - y)J{t,y)dy. 



(5.3) 



Remark. Many equalities in the paper, including (5.3) hold up to a discretization error. To simplify presen- 
tation, we do not mention this in the sequel when discrete sums are approximated by integrals. 
The interaction stress can be rewritten as 



N 



E 



N - 1 



U' 



ipriix - sqj+i ~ {I ~ s)qj)ds. (5.4) 



Next we approximate 



iqj + l-qj)/h^q'it,X) 



{q-^)'{t,q{t,X)) Jit,q{t,X)) 



(5.5) 



This approximation is in fact exact, provided the interpolant is chosen to be piecewise linear. Note also that 
this equation is a special feature of one-dimensional dynamics, where the derivative (deformation gradient) 
can be identified with the Jacobian of the deformation map. In higher dimensions, this no longer holds. 
Inserting this into (5.4), replacing Riemann sum with an integral and changing variable of integration as in 
Sect. 3.2 we obtain the integral approximation 



N - 1 



Jit,y) 



tp^{x-y- 



sh 



ds dy. 



(5.6) 



Equations (5.3 1, (5.6) contain two microscale quantities: J and v. Approximating sums in the definitions of 
the primary averages (2.5), (2.6) by integrals we see that p'' and are obtained by applying the convolution 
operator i?,, to, respectively J and Jv: 



M 
~L 



M 
~L 



(5.7) 
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Recall that Q,, denotes a regularizing approximation to the exact inverse operator R ^ . Applying in 



(5.7) yields integral approximations 

L 



(5.8) 



Inserting these approximation into the integral formulas for the stress yields closed form mesoscopic contin- 
uum equations 



5tp" + a,(p%'') = o, 



where T'L^tI^^. are given by 



T' — 
J (c) - - 



(5.9) 
(5.10) 



(5.11) 



(■ini) 



N - 1 



u' 



M 



Up'']it,y)J Jo 



Ms 



Qr,[p"](t,y) 



ds dy. 



(5.12) 



The choice of the operator is not unique: it depends on the specifics of the regularization method and 
the chosen value of the regularization parameter. For classical regularization schemes such as Landweber, 
Tikhonov, and truncated SVD, this operator is a convolution with the kernel Q{x) that can be described 
explicitly in terms of the singular values and singular vectors of i?,, (see |9] for details). 

In [5D] we studied a special case of (5.8 1 corresponding to the Landweber approximation (3.9 1 with n — 0. 



In that case, called zero-order closure, Q is the identity operator. This means that 



J : 



Vq — v"^ ■ 



(5.13) 



6. Numerical experiments. 

6.1. Lennard-Jones chain. In this example, we simulate a chain of particles interacting with Lennard- 



Jones potential plotted in left panel of Fig. B.l and defined in (B.l ) in the Appendix B. The initial positions 
are equally spaced with 5° — {j — 1/2) ft-, j — 1, . . . ,iV. We consider two different sets of initial velocities 
shown in Fig. 6.1 The left panel contains a meso-scale feature (the larger peak), and a sub-filter scale 
feature (the smaller peak). The right panel shows the same initial velocity but with added noise. The noise 
is a realization of a uniformly distributed random variable. In the sequel, we refer to the first initial condition 
as deterministic, while the second initial condition is called noisy. 
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Fig. 6.1. Left panel: deterministic initial velocity; right panel: noisy initial velocity 



In Fig. |6.2[ we show exact and reconstructed Jacobians. The exact Jacobian in the deterministic case 
differs from the Jacobian in the noisy case, but the reconstructions are similar. The similarity may be due 
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Fig. 6.2. Left panel: reconstruction of the Jacohian J in the deterministic case; right panel: noisy case 

0.04 n 





Fig. 6.3. Left panel: reconstruction of the velocity v in the deterministic case; right panel: noisy case. On both panels 
the exact velocity is shown in red (dark grey) thin solid line, the average velocity in green (light grey) dashed line, and the 
reconstructed velocity in black solid line. 



to the built-in filtering in the deconvolution algorithm. This filtering is rather soft, since the reconstrueted 
Jacobians contain oscillatory artifacts on scales comparable with the micro-scale. The amplitude of the arti- 
facts is under control, so that the relative loo error does not exceed 0.3%. This shows that the reconstruction 
is stable. 

The velocity reconstruction is shown in Fig. |6.3[ The average velocities in the noisy and deterministic 
case are nearly identical. Averaging obliterates sub-lifter scale features, but the deconvolution algorithm 
reco vers these features well. At the same time, the high frequency noise in the right panel is filtered out. In 



Fig. 



6.4 



we show the exact convective stress T^^^ and its closed form approximation T 



computed from 



the equation (5.11). We see that the approximation quality is better in the deterministic case. In the noisy 



case, the main features of the stress are still well recovered, despite a significant difference between the exact 



and reconstructed velocities (Fig. 6.3 1 





Fig. 6.4. Exact convective stress T^'^j (green (grey) dashed line), and the approximation T^^j (black solid line). Left panel: 
deterministic case; right panel: noisy case. 



The interaction stress Tj'.^^^ and its closed form approximation T^^^ (eq. 



5.12)) are shown in Fig. 



6.5 



The approximation quality is about the same in both cases. According to ( |5.12p , the accuracy depends on 
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Fig. 6.5. Exact interaction stress T^^^^^ (green (grey) dashed line), and the approximation T^^j (black solid line). Left 
panel: deterministic case; right panel: noisy case. 

the quality of the reconstruction of the Jacobian. Comparison of Fig. |6.2| and Fig . |6.5| shows that the high 
frequency artifacts and noise in the Jacobian are smoothed out by averaging in ( |5.12[ ). This is important 
since the interaction stress depends non-hnearly on the Jacobian. In contrast to the hnear case, non-hnear 
averaging functionals may exhibit sensitivity to osciUations in the input function (in this case, Jacobian). 
Non-hnearity induces dispersion, and this may result in transfer of the input's high frequency content into 
the low frequency content of the functional. In the present case, this does not happen. We conjecture that 
for generic molecular potentials such as Lennard- Jones, and for a broad class of initial conditions, the stress 
functional has the self-averaging property, meaning that the effect of dispersion is weak compared to the 
filtering effect of convolution. 
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Fig. 6.6. Initial microscopic velocity (black solid curve) defined in fB. 5| ) together with the average velocity v 

(green dashed curve). Left panel: velocities are plotted on the entire domain [0, 1]; right panel: zoom-in of velocities on the 
interval [0.26, 0.34] that contains a sub-filter feature 



6.2. Granular acoustics. 



In this subsection, we test the method on a chain of particles interacting 

The 



B.l 



with a pair potential U (^) defined in the Appendix in (B.3 1 and depicted in the right panel of Fig, 
particles represent the centers of spherical granules and the potential resembles a Hertz potential employed 
in modeling of granular materials. The corresponding force is purely repulsive and has a finite range equal 
to the equilibrium distance between the neighboring particles. 



We solve the system of ODEs (4.3), (4.4) with two different initial conditions and periodic boundary 
conditions. In both examples, the initial positions qj are equally spaced on the interval (0, L) at the equi- 
librium distance h = L/N with L = 1 and = 10,000. The initial velocity for the first example, shown in 
Fig. |6.6| (black curve) and defined in the Appendix (B.4), (B.5|, contains features of different length scales. 
The size of the larger feature of trapezoidal shape is bigger than rjL (mesoscale or filter scale) . The smaller 
feature is a Gaussian with the standard deviation 0.2t]L (sub-filter scale). At this length scale, the feature 



6.6 



is completely obscured by the averaging (green dashed curve) as can be seen in the right panel of Fig, 
where the velocity is zoomed around x = 0.3. 

The system of ODEs (4.3), (4.4) is solved numerically untiH = 2.2-10^^ after which the solution develops 



a shock. First, we test the quality of reconstruction of v and J. Figure [6?7| compares the exact v (red thin solid 
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curve), its reconstructed approximation 



(black thick solid curve) and the average v (green dashed 
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Fig. 6.7. Velocity reconstruction: microscopic velocity v (red thin solid curve), reconstructed velocity ^ (black 

thick solid curve) and average velocity v (green dashed curve). Left panel: on the entire domain; right panel: zoom-in of the 
region that contains a feature due to Gaussian perturbation 



curve) at t = 10 ^. The left panel shows that the reconstruction captures the large scale features (left panel) 
as well as features on the sub-filter scale (right panel). In contrast, the average velocity completely misses 



the sub-filter scale. In Fig. 6.8 we compare the exact microscopic Jacobian J (red thin solid curve) with the 
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Fig. 6.8. Jacobian reconstruction: microscopic Jacobian J (red thin solid curve), reconstructed Jacobian jyQrjlp''] (black 
thick solid curve) and average density -^p (green dashed curve). Left panel: solutions are shown on the entire domain; right 
panel: zoom-in of the region with sub-filter scale features 



reconstructed Jacobian jjQtjIp^] (black thick solid curve) and the average scaled density -^p (green dashed 



curve). Similar to velocity reconstruction, Fig. 6.8 indicates that the reconstructed Jacobian is much closer 
to the exact Jacobian than the average density. 

stresses are approximated by T'A-i and 



Next we examine how well convective T\ and interaction TP. 



6.9 



indicates that the exact convective stress and 



T(;„t) defined in (5.11), (5.12). The left panel of Fig. 

its approximation are alniost indistinguishable. The right panel of Fig. |6.9| shows good agreement between 
the exact interaction stress and its approximation. For comparison, we also plot an approximation using a 



zero-order closure from shown for convenience in (5.13) that fails to capture sub- filter scale features in 
both stresses. 



The L 



-error in approximation of T^'^^ by T^^^ is between 1.5% and 10% during the simulation time. The 



error in approximation of is smaller and varies from 1.5% to 8%. Preliminary computational studies 

indicate that the error decreases as N increases. 

The initial velocity in the second example is the sum of u^"* 



defined in (B.4) and a sine function with 



period 0.012, added on the interval [0, 0.6] (see Fig. 6.10 and formu la in (B.6 1 in the Appendix). Simulations 
were also done until t = 2.2 • 10~^. The right panel of Fig. 6.10 shows that at t — the average velocity 
does not contain oscillations present in the microscopic velocity. Fig. |6.11| presents graphs of velocity and 
Jacobian at a representative moment of time, t — 10^'^. The reconstructed velocity and Jacobian contain 
main features of their microscopic counterparts while the averages do not. 



Fig. 6.12 depicts the stress. Both convective and interaction exact stresses have sharp transition regions 
near x — 0.05 and x = 0.6. Our closed form approximation qualitatively captures these features while a 
zero-order closure approximation is nearly zero on the entire interval. The error in approximation of the 
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Fig. 6.9. Left panel shows convective stresses: exact T^^^ 
(black thick solid curve), and an approximation (green dash 
stress T^. , and its corresponding approximations 
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Fig. 6.10. Initial velocity with sine perturbation: micros: 
dashed curve). Left panel: velocity is shown on the entire dorr. 
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(red thin solid curve), its approximation via reconstruction T^^^ 
curve) using zero-order closure (KT^); right panel: interaction 
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ic velocity (black solid curve) and average velocity v (green 
right panel: zoom-in of the velocity on the interval [0.5, 0.65] 



convective stress fluctuates at early times (until t — i ■ 10^^) from 15% to 70% and then settles around 
35 — 40%. The error in the approximation of the interaction stress behaves similarly at early times, then 
decreases to 10% and levels off. The error in using the zero-order closure is much higher: 75 — 100% for 
the convective stress during the entire simulation time and around 100% for the interaction stress until 
t = 2 ■ 10"'^, then it drops to 10 — 15% at t = 7 ■ 10^^, after which the error is about the same as using the 
reconstruction. 

7. Conclusions. We propose a method for deriving closed form mesoscale continuum models of large 
particle systems. The closure construction is based on the following. Non- linear meso-scale averages can be 
rewritten as linear convolutions of the window function and appropriate micro-scale dynamical functions. One 
such function of particular importance is the inverse Jacobian of the micro-scale flow map associated with a 
position interpolant. Using the theory of ill-posed problems, we produce stable deconvolution approximations 
of particle positions and velocities in terms of the average density and average momentum. Closure is achieved 
by inserting these approximations into the equations for fluxes instead of the actual positions and velocities. 



The resulting constitutive equations (5.111, (5.12) are non-local in space and non-linear. 

In the simplest version of the method, the micro-scale quantities are approximated by their averages. We 
studied this approximation in the earlier paper |20j . The results presented there indicate that the simplest 
approximation works well for systems characterized by (i) small fluctuations of the initial velocity; and (ii) 
nearly isothermal dynamics. In this article we consider more general initial conditions that contain prominent 
small scale peaks, significant noise, or high frequency periodic oscillations. Averaging obscures these features 
to such an extent that the approximation from |20| becomes unsatisfactory. Here we were able to recover 
such details using non-iterative regularization methods. 

We tested the method numerically on two models of FPU-chains: the classical Lennard-Jones chain, 
and the granular acoustics model considered earlier in |20j . but with more general initial conditions. The 
ODEs were solved by the velocity Verlet method, and the obtained particle positions and velocities were 
used to calculate the average density, linear momentum, and the exact stress. Then we used regularized 
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Fig. 6.11. Reconstruction of velocity v (left panel); reconstruction of the Jacobian (let panel) at t = 10~^. The functions 
are plotted on the interval [0.5, 0.65] to show details. 
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Fig. 6.12. Left panel: convective stress; right panel: interaction stress 

dcconvolution to generate the approximate (reconstructed) Jacobian and velocity. The resuhing closed form 
approximation of the stress agreed very well with its exact counterpart. 
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Appendix A. Window function. 

In this paper, we use the following function ip: 



m = { 



a+b ' 
_A+b 



I 0, 



if iei<a, 

if a < ^ <b, 

b2, if ~b<£,<-a, 

if m>b 



(A.l) 



with L = 1, a = L/2, b = 3L/2. It can be directly checked that "0(0^^^ = /_b i^iO'^^ ~ f • The function 
i/) is plotted in Fig. |A.1| 

Appendix B. Potentials and initial conditions. 

In Section [6j we test the method using two potentials: the first is the Lennard- Jones potential, the 
second potential is similar to the Hertz potential employed in modeling of granular materials. 



B.l. Lennard- Jones. The potential is plotted in the left panel of Fig. B.l and defined by 

12 / \ 6' 



(B.l) 



where e = 0.25 defines the depth of the potential well, a is the finite distance at which the potential is zero, 
^ is the distance between particles. The potential is at a minimum when ^ = h — 2^/^cr, which determines 
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Fig. A.l. The function 

the choice of a. The force corresponding to the Lennard- Jones potential is repulsive for distances smaller 
than h and attractive for distances greater than h. 
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Fig. B.l. Left panel: Lennard- Jones potential; right panel: Hertz potential 

The initial particles positions in all numerical test problems are equally spaced with step h and defined 
by q'j — {j — l/2)h, j — 1, . . . , N. The initial velocity used in the deterministic case is 



^^(9°) - /('Zj) + A - 0.7) , j = l,...,N (B.2) 

where 



/(0 = 



and 



5b(?-i) (1-0 if i<e<i, 
otherwise. 



otherwise. 



The initial velocity for the noisy case is the sum of vq from (B.2) and a uniformly distributed random 
variable with mean zero and maximum amplitude 10""^. 

B.2. Granular acoustics. The potential is defined as 

\ 0, a 

where p > 1, x^, — L, and Cr is material stiffness. The potential is plotted in the right panel of Fig. |B.1| 
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The initial velocity in the first example of Subsection 6.2 is given by = ^base _j_ ^i,pert^ where v^"'"^ is 
a piecewise cubic continues function and z;^ is a Gaussian: 



base / ^0 



V 



l^pert / 



0, 




if 


< g," < ii, 




-xi)(g^"-Li)2, 


if 


Li < <z° < L2 


d2, 




if 


L2 < q° < L3 






if 


L3 < q° < U 


0, 




if 


Li < q° < L, 



(g ) = ai exp 



iq^~q*y 

2cr2 



J = 1, . 



J = 1, . 



.,iV, 



(B.4) 



(B.5) 



Here Li = 0.2L, L2 = 0.4L, L3 0.7L, L4 ^ 0.9L, xi = (3L2-ii)/2, X2 = (3L3-i4)/2, ^2 = 0.3, 



di 



-2d2/[L2 - Lif, ds = -2d2/(L3 - Lif, = 0.1, q* = 0.3L. 



The initial velocity in the second example is 



V 



base 



+ v^^pert^ where u^"^*^ is as in (B.4 1 but with 



Li = O.IL, L2 = 0.2L, L3 = 0.3L, L4 = 0.6i and ti^^P^^'* is a sine function on the interval [0, L4]: 

' a2sin('^), if ^<q]<U, J = 1, . 
0, otherwise 



V 



2, pert / 



,N, 



(B.6) 



with 02 — 5 and fc = 50. The sine perturbation has period 0.012. 



REFERENCES 

[1] N. A. Adams and S. Stolz, A subgrid-scale deconvolution approach for shock capturing, J. Comp. Phys., 178 (2002), 
pp. 391-426. 

[2] L. C. Berselli, T. Iliescu, and W. J. Layton, Mathematics of Large Eddy Simulation of Turbulent Flows, Springer, 
New York, 2006. 

[3] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer Academic, Dordrecht, 1996. 
[4] V. Fridman, a method of successive approximations for Fredholm integral equations of the first kind, Uspekhi Mat. Nauk, 

11 (1956), pp. 233-234 (in Russian). 
[5] C. W. GrOETSCH, The Theory of Tikhonov Regularization for Fredholm Equation of the First Kind, Pitman, Boston, 

1984. 

[6] P. Ch. Hansen, Rank- Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion, SIAM, 1987. 
[7] R. ,J. Hardy, Formulas for determining local properties in molecular- dynamics simulations: shock waves, J. Chem. Phys., 
76 (1982), pp. 622-628. 

[8] J. H. Irving and J. G. Kirkwood, The statistical theory of transport processes IV. The equations of hydrodynamics, J. 

Chem. Phys., 18 (1950), pp. 817-829. 
[9] A. KiRSCH, An Introduction to the Mathematical Theory of Inverse Problems, Springer, New York, 1996. 
[10] L. Landweber, An iteration formula for Fredholm integral equations of the first kind. Am. J. Math., 73 (1951), pp. 615— 

624. 

[11] W. Layton, Bounds on helicity and dissipation rates of approximate deconvolution models of turbulence, SIAM J. Math. 
Anal., 39 (2007), pp. 916-931. 

[12] W. Layton and O. Lewandowski, Residual stress of approximate deconvolution models of turbulence, J. Turbulence, 7 
(2006), pp. 1-21. 

[13] R. B. LehouCQ and M. P. Sears, The statistical mechanical foundation of the peridynamic nonlocal continuum theory: 

energy and momentum conservation laws, Phys. Rev. E, (to appear). 
[14] V. A. MOROZOV, Methods for Solving Incorrectly Posed Problems, Springer, New York, 1984. 
[15] A. I. Murdoch, A critique of atomistic definitions of the stress tensor, J. Elasticity, 88 (2007), pp. 113-140. 
[16] A. I. Murdoch and D. Bedeaux, Continuum equations of balance via weighted averages of microscopic quantities, Proc. 

Royal Soc. London A, 445 (1994), pp. 157-179. 
[17] , A microscopic perspective on the physical foundations of continuum mechanics - Part I: macroscopic states, 

reproducibility, and macroscopic statistics, at prescribed scales of length and time. Int. J. Engng Sci., 34 (1996), 

pp. 1111-1129. 

[18] , A microscopic perspective on the physical foundations of continuum mechanics - Part II: a projection operator 

approach to the separation of reversible and irreversible contributions to macroscopic behaviour. Int. J. Engng Sci., 
35 (1997), pp. 921-949. 

[19] W. Noll, Der herleitung der grundgleichungen der thermomechanik der kontinua aus der statistischen mechanik, J. 
Ration. Mech. Anal., 4 (1955), pp. 627-646. 



19 



[20] A. Pancmenko, L. L. Barannyk, and R. P. Gilbert, Closure method for spatially averaged dynamics of particle chains, 

Nonlinear Anal. Real World Appl., 12 (2011), pp. 1681-1697. 
[21] G. A. Pavliotis and A. M. Stuart, Multiscale methods. Averaging and Homogenization, Springer, 2008. 
[22] S. Silling and R. B. Lehoucq, Peridynamic theory of solid mechanics, Advances in Applied Meehanics, 44 (2010), 

pp. 7.3 168. 

[23] A. Tartakovsky, A. Panchenko, and K. Ferrjh, Dimension reduction method for ODE fluid models, J. Comp. Phys., 

(to appear), p. Preprint available at littp://www.niatli.wsu.cdu/inath/faculty/panchenko/welcome.php. 
[24] A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, Wiley, New York, 1987. 



20 



